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Preface 

Computer  modeling  is  extensively  applied  in  all  areas  of  engineering  and  science.  Especially 
for  the  design  of  economical  and  efficient  aerospace  structures  computer  modeling  is  a 
valuable  tool.  For  flow  around  aerospace  structures  large  computer  memory  and  time  is 
required  due  to  complex  geometry.  If  accurate  solution  is  obtained  with  less  number  of  grid 
point  that  is  the  ideal  situation.  To  achieve  this  objective  adaptive  Finite  Element  Method 
(FEM)  is  tried  here. 

The  adaptive  FEM  technique  was  proposed  in  early  1980s  to  reduce  the  discretized 
errors  resulted  from  a  FEM  mesh.  In  this  technique,  the  errors  are  first  calculated 
to  assess  the  accuracy  of  the  solution.  If  the  errors  are  large,  the  finite  element 
model  is  then  refined  through  redistributing  the  nodes  (called  R-version  adaptive 
FEM),  or  reducing  the  size  of  elements  (called  H-version  adaptive  FEM),  or 
increasing  the  order  of  the  interpolation  functions  (called  P-version  adaptive  FEM) 
[13],  or  using  the  combination  ways.  The  new  model  is  then  re-analyzed  and  the 
errors  in  the  new  model  are  recalculated.  The  procedure  is  continued  until  the 
calculated  errors  fall  below  the  specified  permissible  values.  Hence,  the  adaptivity 
means  that  the  FEM  model  refinement  is  based  on  the  error  distribution. 

For  standard  problem  with  regular  geometry  several  work  has  been  done  using 
adaptive  FEM  technique.  However,  there  seems  to  be  far  less  work  being  done 
directly  for  complex  geometries.  There  are  several  challenges  like  selection  of  proper 
solvers  when  the  condition  number  is  very  high  and  the  advantage  of  using  p- 
adaptive  and  h-adaptive.  These  issues  are  addressed  in  this  research. 

We  would  like  to  express  our  special  thanks  to  the  AFOSR,  DOD  for  their 
sponsorship  of  this  project.  Especially  the  encouragement  provided  by  Dr.  Len 
Sakell,  AFOSR  in  the  initial  stages  is  acknowledged. 

R.  Panneer  Selvam  &  Zu-Qing  Qu 
Fayetteville,  Arkansas,  May  2004 
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Abstract 


Predicting  wind  induced  aerodynamic  forces  on  aerospace  structures  is  usually 
computationally  expensive  even  for  powerful  computational  facilities.  For  these 
unsteady  computations,  adaptive  finite  element  technique  may  reduce  the 
computer  time  and  storage  while  keeping  a  desired  accuracy.  In  this  report,  the  p- 
version  adaptive  finite  element  method  is  implemented  into  a  standard  benchmark 
problem,  the  computational  flow  around  a  circular  cylinder,  to  compute 
aerod5Tiamic  forces.  The  error  distribution  for  velocity  is  first  estimated.  Then,  the 
polynomial  order  of  the  interpolation  function  is  changed  continuously  in 
accordance  with  the  changing  of  the  error  distribution.  The  second  through  fifth 
orders  of  polynomials  are  considered  for  the  velocity  in  the  adaptive  method.  One 
order  less  of  polynomial  is  used  for  pressure.  The  benchmark  problem  of  the  flow 
around  a  circular  cylinder  with  Reynolds  number  of  1000  is  considered  to  study  the 
performance.  The  effects  of  the  highest  order  of  polynomials,  error  tolerance,  and 
size  of  the  element  on  the  accuracy  of  the  drag  and  lift  coefficients  are  surveyed 
using  this  flow  simulation.  The  results  show  that  the  accuracy  of  the  velocity  close 
to  the  cylinder  affects  the  drag  and  lift  coefficients  and  the  error,  20%  for  example, 
far  away  from  the  cylinder  does  not  have  much  effect  on  the  accuracy  of  these 
coefficients. 


1.  Introduction 

r 

Prediction  of  wind-induced  aerodynamic  forces  on  aerospace  structures  has  been 
mainly  made  by  experimental  methods  involving  the  wind  tunnel  technique  [1,2]. 
However,  with  the  advent  of  supercomputers  the  tendency  is  gradually  changing. 
Powerful  computationed  facilities  make  it  possible  to  deal  with  these  physical 
quantities  numerically.  The  Large  Eddy  Simulation  (LES)  based  on  the 
Smagorinsky’s  eddy  viscosity  model  used  by  Murakami  and  his  research  group  for 
wind  engineering  [3]  is  used  in  the  modeling  of  drag  crisis  around  a  circular 
cylinder.  Kakuda  and  Toskada  [4],  Kondo  [5]  and  Tamura  et  al  [6,7]  used  the  no- 


4 


Adaptive  Navier-Stokes  Flow  Solver  for  Aerospace  Structures 


turbulence  model  to  study  the  flow  around  a  circular  cylinder.  Instead  of  using 
explicit  methods  [3-7],  an  implicit  solution  procedure  to  solve  the  Navier-Stokes  (NS) 
equations  using  LES  and  finite  element  method  was  proposed  by  Selvam  [8]  to 
study  flow  around  a  cylinder.  As  we  know,  to  solve  the  flow  around  a  bluff  structure 
is  highly  computationally  expensive.  For  these  unsteady  computations,  adaptive 
finite  element  method  may  reduce  the  computer  time  and  storage. 

The  adaptive  FEM  technique  was  proposed  in  early  1980s  to  reduce  the  discretized 
errors  resulted  from  a  FEM  mesh.  In  this  technique,  the  errors  are  first  calculated 
to  assess  the  accuracy  of  the  solution.  If  the  errors  are  large,  the  finite  element 
model  is  then  refined  through  redistributing  the  nodes  (called  R-version  adaptive 
FEM)  [9,10],  or  reducing  the  size  of  elements  (called  H-version  adaptive  FEM) 
[11,12],  or  increasing  the  order  of  the  interpolation  functions  (called  P-version 
adaptive  FEM)  [13],  or  using  the  combination  ways.  The  new  model  is  then  re¬ 
analyzed  and  the  errors  in  the  new  model  are  recalculated.  The  procedure  is 
continued  until  the  calculated  errors  fall  below  the  specified  permissible  values. 
Hence,  the  adaptivity  means  that  the  FEM  model  refinement  is  based  on  the  error 
distribution. 

The  adaptive  FEM  technique  has  been  discussed  in  detail  during  the  past  two 
decades.  The  results  of  such  work  are  fruitful  [14].  However,  there  seems  to  be  far 
less  work  being  done  directly  in  the  wind  engineering  which  needs  much  more 
complicated  procedures. 

Choi  and  Yu  [15,16]  investigated  the  h-refinement  for  flows  over  a  square  cylinder. 
They  used  the  penalty-function  formulation  to  solve  the  NS  equations.  This 
procedure  is  not  used  commonly.  Selvam  [17,  18]  applied  the  mesh  enrichment 
technique  (h-refinement)  and  p-refinement  techniques  to  flows  over  a  circular 
cylinder.  He  used  primitive  variable  form  to  solve  the  NS  equations. 

In  this  report,  the  finite  element  modeling  of  the  benchmark  problem,  the  flow 
around  a  circular  cylinder,  is  described  in  Section  2.  Three- step  advancement 
scheme  for  solving  the  N-S  equations  by  LES  is  discussed.  The  error  estimation 
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based  on  the  velocity  is  proposed  in  Section  3.  The  details  on  the  implementation  of 
the  p-adaptive  technique  are  also  discussed  in  this  section.  In  Section  4,  the 
simulation  of  a  flow  over  a  circular  cylinder  for  Reynolds  number  of  1000  is 
considered.  The  effects  of  the  highest  order  of  polynomials,  error  tolerance,  and  size 
of  the  element  on  the  accuracy  of  the  drag  and  lift  coefficients  are  surveyed  using 
this  flow  simulation.  Some  useful  conclusions  will  be  given  in  that  section.  The 
computed  drag  and  lift  coefficients  are  also  compared  with  the  experimental  and 
other  computational  results. 


2.  Finite  Element  Modeling  of  the  Flow 

In  the  following  discussion  Reynolds  number  ,  drag  coefficient  Q ,  lift  coefficient 

C, ,  and  Strouhal  number  S,  are  defined  as 

'R.^VDlv 

C,=Fj{0.5py^D) 

‘  C,=Fj(0.5pV’D) 

S,=DfTV 

where  D  is  the  diameter  of  the  cylinder,  V  is  the  reference  velocity,  v  is  the 
kinematic  viscosity,  and  Fy  are  the  drag  and  lift  force,  T  is  the  period  of  oscillation 

of  the  Hft  forces  and  p  is  the  density. 

The  flow  around  a  cylinder  is  represented  by  using  the  Navier-Stokes  (NS) 
equations.  The  two  and  three-dimensional  equations  for  ah  incompressible  fluid  in 
general  tensor  notation  are  as  follows: 

Continuity  Equation:  =  0  (2) 

Momentum  Equation:  C/, ,  +  Up^j  =  -{pf  yo)  ,•  +  [v(c/,-  ^  (3) 

where  t/,.  and  p  are  the  velocity  and  pressure  respectively,  p  is  the  fluid  density.  A 

comma  represents  one  differentiation;  t  represent  time,  i  =1,2  and  3  mean  variables 
in  the  x,  y  and  z  directions.  To  implement  higher  order  approximation  of  the 
convection  term  the  following  expression  is  used  instead  of  Equation  (3) 
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t/,, + up,,  -  e[upp,,,),  /2 = -(pIpI + [v(u,j + u,,  )J 


(4) 


Depending  upon  the  values  of  6  ,  different  procedures  can  be  implemented.  For 
balance  tensor  diffusivity  (BTD)  scheme  0  =  is  used;  where  5t  is  the  time  step  size 
used  in  the  integration.  For  streamline  upwind  procedure  suggested  by  Brooks  and 
Hughes  [19],  0  is  considered  as 


max 

[dx  dy  dz ) 

Here  dx,  dy  ,  and  dz  are  the  control  volume  length  in  the  x ,  y  ,  and  z  directions;  C/, , 
U2,  and  f/3  are  the  velocities  in  the  three  directions.  In  this  computation  Q  =bt  is 
used.  This  has  less  numerical  diffusion  as  compared  with  benchmark  problems  [20]. 


1 


fc/,  V,  \ 

U,] 

.dx’dy’ 

dz  j 

The  NS  equations  are  solved  by  using  an  implicit  method,  which  is  similar  to 
Selv£un  [8],  to  eliminate  the  numerical  stability  restrictions.  The  three-step 
advancement  scheme  for  Equations  (2)  and  (4)  is  as  follows: 


Step  1:  Solve  for  [/,.  from  Equation  (3).  The  diffusion  and  the  higher  order 

convection  terms  are  considered  implicitly  to  be  in  the  current  time  and  the  first 
order  convection  terms  are  considered  explicitly  from  the  previous  time  step.  The 
pressure  is  considered  in  the  right  hand  side  of  the  equation.  This  set  of  equations 
leads  to  a  S5mimetric  matrix  and  the  preconditioned  conjugate  gradient  (PCG) 
procedure  is  used  to  solve  it. 

Step  2:  Solve  for  pressure  correction  from 

(^X  =  u„ia. 


Step  3:  Correct  the  velocity  for  incompressibility: 

u,^u,-a(sp,) 

where  C/,  is  not  specified  and  update  the  pressure  p  =  p  +  5p  . 
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Suppose  the  element  variables  and  P  are  discretized  using  the  shape  functions 


Ny  and  N ^  ,  i.e. 

U]  =  NyU,,  P'=N^p  (6) 

where  Ny  and  are  1  x  vectors;  «,  and  p  are  x  1  vector;  is  the  number  of 

nodes  for  each  element.  Here,  the  velocity  and  pressure  are  approximated  by  using 
unequal  order  interpolation  polynomial  because  if  an  equal  order  is  used,  the 
solution  starts  to  diverge  due  to  the  violation  of  Babusk  and  Brezzi  condition. 
Introducing  Equation  (6)  into  Equation  (2)  and  using  the  regular  Galerkm 
procedure  the  following  matrix  expressions  are  obtained: 

WK,. ,  +  (</  +  V,  +  Vj )«,  +qp  =  f  (7) 


where 


In  equation  (7),  m  ,  d  ,  v, ,  Vj,,  and  q  are  the  mass,  dififusion,  convection  due  to  the 
first  and  second  order  and  pressure  matrices  of  size  n^xn^.  f  is  the  n^xX  vector 

which  considers  the  given  Neumann  boundary  conditions.  After  assembling  all  the 
element  matrices,  the  equations  of  the  finite  element  model  becomes 

MU,,+{D-¥V,+V^)U,^QP  =  F  (9) 

Equation  (9)  can  be  solved  by  using  the  aforementioned  implicit  method.  The 
detailed  solution  procedure  is  similar  to  Selvam  [8]. 


3.  P-adaotive  Finite  Element  Technique 

The  error  estimate  is  one  of  the  most  important  steps  in  the  adaptive  technique.  It 
gives  the  error  distribution  in  the  present  finite  element  grid,  which  can  be  used  as 
indicator  to  refine  the  grid.  The  discretization  errors  of  a  finite  element  solution  can 
be  estimated  by  implicit  or  explicit  method  [21].  In  structural  engineering,  stress 
recovery  techniques  are  usually  used  to  estimate  the  errors  because  they  are  much 
easier  to  implement  in  the  programming.  As  the  exact  solution  is  generally  not 
known,  most  of  the  approaches  are  concerned  with  posteriori  error  estimates. 
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In  the  present  research,  the  error  estimator  based  on  the  velocity  is  considered.  The 
error  is  computed  by  considering  the  higher  order  solution  to  be  exact  and  the  lower 
order  to  be  approximate.  The  difference  between  the  two  is  defined  as  the  error  in  the 
present  grid,  that  is 


(10) 


where  v^.,  and  are  the  velocity  obtained  from  the  finite  element  model  using  p-\ 
and  p  order  interpolation  functions,  respectively.  The  order  of  polynomial  to  start 

with  is  quadratic  (p=2)  for  velocity.  Hence  it  is  easy  to  compute  firom  quadratic  to 
linear,  third  order  to  quadratic,  and  fourth  order  to  third  order. 


After  the  errors  have  been  estimated  and  been  found  to  be  big,  the  next  step  is  to 
refine  the  finite  element  model  so  as  to  reduce  the  errors.  The  accuracy  of  a  finite 
element  solution  depends  upon  the  shape  and  size  of  the  elements  and  the  order  of 
the  interpolation  functions.  Consequently,  there  are  three  methods  to  refine  the 
finite  element  model.  In  this  paper,  p-version  is  implemented  to  the  wind 
engineering. 

P-refinement  increases  the  order  of  the  interpolation  functions  while  keeping  the 
mesh  unchanged.  Higher  order  elements  generally  provide  a  better  description  of 
the  domain  geometrically.  They  are  particularly  useful  in  regions  where  use  of  lower 
order  elements  would  result  in  a  mesh  with  poor  aspect  ratios  in  those  elements 
[14].  From  the  point  of  view  of  solution  accuracy,  higher  order  elements  are  usually 
more  accurate  than  the  lower  order  elements. 


The  Hierarchical  functions  are  applied  here  to  increase  of  the  interpolation  function 
order  because  the  method  requires  only  the  computation  of  the  coefficients  in  rows 
and  columns  associated  with  the  new  enriched  degrees  of  freedom,  which  together 
with  the  previously  computed  element  matrices  form  a  new  stiffness  matrix  [22,23]. 

For  one-dimensional  elements,  the  linear  shape  functions  are  usually  defined  as 
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^ofe)  = 


1-^ 

2  ’ 


(11) 


where  ^  (-l<^<l)isa  non-dimensional  coordinate.  There  are  many  ways  to 
construct  the  higher  order  shape  function  [23].  The  following  shape  functions  are 
used  in  this  report: 


J^s(0  = 


J^*-l  5  even 
1^*  -(^  s  odd 


(12) 


where  s  (>2)  is  the  degree  of  the  introduced  polynomial.  The  corresponding 
displacement  function  can  be  obtained  as 


1^0 

Using  these  one-dimensional  formulae,  it  is  easy  to  derive  two-dimensional  shape 
functions  [18].  These  shape  functions  can  be  directly  used  in  Equations  (8)  to 
compute  the  element  matrices. 


In  the  present  work,  four-node  quadrilateral  element  is  used  to  describe  the 
geometry  of  the  element.  The  polynomial  order  considered  for  velocity  are  from  2  to 
5,  which  means  in  each  direction  interpolation  function  of  the  order  of  2  through  5 
are  considered.  The  1st  through  4th  orders  of  polynomial,  which  are  one  order  less 
than  that  for  velocity  respectively,  are  considered  for  the  pressure.  The  size  of  the 
element  stiffness  matrix  for  velocity  varied  from  9x9  to  36x36.  The  matrices  are 
numerical  integrated.  The  integration  points  considered  at  this  time  are  3x3 
through  36x36. 


4.  Simulation  of  the  Flow  around  a  Circular  Cylinder 

The  computational  region  of  the  wind  around  a  circular  cylinder  is  shown  in  Figure 
1.  Its  length  in  the  a:  and  y  directions  are  non-dimensionalized  with  respect  to  the 
diameter  of  the  cylinder.  The  inlet  velocities  in  these  two  directions  are  considered 
to  be  one  and  zero,  respectively.  On  the  top  and  bottom  sides,  the  vertical  velocities 
and  the  normal  derivative  of  the  velocities  are  set  to  be  zero.  The  velocities  are  also 
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considered  to  be  zero  (no  slip)  on  the  surface  of  the  cylinder.  The  flow  around  the 
circular  cylinder  for  Reynolds  number  of  1000  will  be  considered  in  the  following. 


1 

2 

Yf 

X 

. .  41  * 

Figure  1.  Computational  region  of  the  wind  around  the  cylinder 

4.1  Coarse  Grid 

Two  coarse  grids  shown  in  Figure  2  are  considered  at  first.  Both  grids  have  414 
(24x15  +  9x6)  elements  and  447  nodes.  The  minimum  spaces  in  the  radial  direction 
around  the  cylinder  are,  respectively,  0.098  and  0.020  in  Grid  A  and  Grid  B.  The 
highest  orders  of  the  interpolation  polynomials  considered  are  P  =  2  and  P  =  5. 
Here  and  in  the  following,  P  denotes  the  highest  pol5momial  order  used  in  the 
adaptive  refinement.  P=2,  for  example,  means  that  no  refinement  is  applied  even 
though  the  error  is  larger  than  the  prescribed  one.  The  flow  is  run  60  seconds.  After 
3  seconds,  the  error  distribution  will  be  evaluated  every  0.2  second.  Tlie  finite 
element  model  will  be  refined  if  the  estimated  error  is  higher  than  prescribed  which 
is  set  as  10%  at  this  time. 


Figure  2.  Two  coarse  grids:  (a)  Grid-A;  (b)  Grid-B 
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(c)  (d) 

Figure  3.  Distribution  of  the  polynomial  order:  (a)  Grid-A;  (b)  Grid-B; 
(c)  Zoom  of  Grid-A;  (d)  Zoom  of  Grid-B 


The  distributions  of  the  polynomial  order  at  the  end  of  60-second  are  plotted  in 
Figures  3(a)  and  3(b).  The  corresponding  error  distributions  are  shown  in  Figure  4. 
The  estimated  errors  are,  respectively,  24%  and  19%  for  the  two  cases.  For  most  of 
the  elements  around  the  cylinder  and  on  its  right  side,  their  orders  have  been 
increased  by  3,  4,  and  5  respectively.  The  5th  order  of  polynomial  are  only  located 
on  the  vorticity  shade  area.  For  clarity  purpose,  the  areas  around  the  cylinder  are 
zoomed  in  Figures  3(c)  and  3(d)  respectively.  Since  the  size  of  the  element  in  the 
radial  direction  used  in  Grid-A  is  much  bigger  than  the  Grid-B,  the  orders  of  many 
elements  around  the  cylinder  in  the  former  case  cire  increased  by  5,  while  they  are  2 
through  4  in  the  latter  case.  As  we  know,  if  the  orders  of  all  the  element  are 
increased  from  2  to  5,  the  number  of  the  unknowns  will  be  increased  by  4  times, 
while  they  are  2.95  and  2.64  times  for  the  Grid-A  and  Grid-B  respectively. 
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(c)  (d) 

Figure  4.  Distribution  of  the  error;  (a)  Grid- A;  (b)  Grid-B; 
(c)  Zoom  of  Grid-A;  (d)  Zoom  of  Grid-B 


The  drag  and  lift  force  coefficients  for  both  cases  are  shown  in  Figures  5  and  6 
respectively.  The  results  obtained  from  both  P  =  2  and  P  =  5  are  included.  The 
difference  of  the  coefficients  is  significant  especially  for  the  Grid-A.  The  average, 
amplitude  of  the  drag  coefficients  and  the  period  of  lift  coefficient  for  the  two  cases 
are  listed  in  Table  1.  The  averaged  drag  coefficients,  for  example,  in  Figure  5  (a)  are 
0.872  and  1.298  for  P  =  2  and  P  =  5  respectively.  The  former  is  much  smaller  than 
the  reasonable  value  1.42  which  will  be  provided  later.  This  means  Grid-A  is 
unreasonable  to  be  used  to  perform  the  force  coefficient  analysis.  After  the  p- 
adaptive  technique  is  implemented,  the  averaged  drag  coefficient  increases  a  lot.  Its 
error  reduces  from  38.6%  to  8.5%.  The  amplitudes  of  the  drag  coefficient  have 
similar  phenomenon. 
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(a)  (b) 

Figure  5.  Drag  coefficient  for  both  cases:  (a)  Grid-A;  (b)  Grid-B 


(a)  (b) 

Figure  6.  Lift  coefficient  for  both  cases:  (a)  Grid-A;  (b)  Grid-B 


Table  1 .  Drag  coefficients  and  period  for  Grid-A  and  Grid-B 


Grid 

P 

Average  (c^) 

Amplitude  (c^) 

Period  ( ) 

Grid- 

2 

0.872 

0.074 

4.35 

A 

5 

1.298 

0.146 

5.00 

Grid- 

2 

1.208 

0.100 

5.00 

B 

5 

1.329 

0.165 

5.00 

Since  the  minimum  space  around  the  cylinder  in  the  radial  direction  is  reduced 
from  0.098  in  Grid-A  to  0.020  in  Grid-B,  the  computed  drag  coefficient  increases  to 
1.208  for  P  =  2.  After  the  refinement  is  considered,  the  coefficient  becomes  1.329. 
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The  difference  between  P  =  2  and  P  —  5  becomes  9.1%  while  it  is  32.8%  in  Grid-A. 
One  of  the  reasons  may  be  that  the  orders  of  the  elements  around  the  cylinder  for 
Grid-B  are  not  increased  as  much  as  that  in  Grid-A.  Totally,  the  estimated  errors 
are  still  a  little  high,  around  20%,  even  though  P  =  5  is  considered  in  Grid-B. 
Therefore,  the  value  of  the  drag  coefficient  1.329  is  a  little  smaller  than  1.42. 
Consequently,  more  refine  grid  is  necessary  for  accuracy  purpose.  The  sizes  of  the 
element  around  the  cylinder  will  be  reduced  both  in  the  radical  direction  and  its 
tangential  direction. 

4.2  Refined  Grid 

Three  finite  element  grids  shown  in  Figure  7  are  considered.  The  properties  of  them 
are  listed  in  Table  2.  The  flow  is  also  run  60  seconds.  After  3  seconds,  the  finite 
element  model  will  be  refined  every  one  second. 


Figure  7.  Three  finite  element  grids;  (a)  Grid-I;  (b)  Grid-II;  (c)  Grid-Ill 
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Table  2.  Properties  of  the  three  finite  element  grids 


Grid  No. 

No.  of  Elements 

No.  of  Nodes 

Eqs.  ofV/P 

Min.  Space 

Grid-I 

1995(60x31+15x9) 

2064 

8118/2064 

0.020 

Grid-II 

3195(75x39+30x9) 

3279 

12948/3279 

0.015 

Grid-Ill 

7335(90x39+45x45) 

7470 

29610/7470 

0.012 

Effects  of  the  Hiehest  Polynomial  Order  P 

The  drag  and  lift  force  coefficients  computed  using  Grid-I  are  plotted  in  Figures  8 
and  9.  The  drag  coefficients  between  40-second  and  60-second  are  zoomed  in  the 
plot  7(b)  for  clear  purpose.  The  prescribed  error  limitation  is  set  as  10%  at  first. 
After  the  simulation  runs  40  seconds,  both  the  drag  and  the  lift  force  coefficients 
become  stable.  It  cein  be  seen  from  Figures  8  and  9  that  there  is  no  much  difference 
among  the  amplitudes  of  the  drag  and  lift  coefficients  for  different  Ps.  For  accuracy 
purpose,  only  the  last  four  periods  of  drag  plots  and  three  periods  of  lift  plots  are 
considered  for  comparison  in  the  following. 


(a)  (b) 


Figure  8.  Drag  coefficients  for  P=2,  3,  4,  and  5  with  Grid-I:  (a)  Full  plot;  (b)  Locally 

zoomed  plot 
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Figure  9.  Lift  coefficients  for  P=2,  3,  4,  and  5  with  Grid-I 

The  averages,  amplitudes,  and  periods  of  the  drag  and  lift  force  coefficients 
computed  from  the  Grid-1  are  listed  in  Table  3.  When  the  polynomial  orders  of  all 
the  elements  are  two,  the  average,  amplitude,  and  period  of  the  drag  coefficient  are 
1.386,  0.195,  and  2.389  respectively.  When  the  highest  poljmomial  order  P  is  set  as 
3,  the  average  and  period  decrease  to  1.370  and  2.369  while  the  amplitude  of  the 
drag  coefficient  increase  a  little.  After  that,  even  though  the  highest  polynomial  P 
has  been  increased  to  4  suid  5,  the  average,  amplitude,  and  period  change  very 
slightly.  As  for  the  lift  force  coefficient,  it  changes  a  little  with  the  increase  of  the 
highest  pol5momial  order  P.  When  the  prescribed  error  is  5%,  these  results  change  a 
lot  from  P=2  to  P=3.  Then,  there  is  no  much  difference  for  the  P=3,  4,  and  5. 


Table  3.  Drag  and  lift  coefficients  for  Grid-I 


Prescribe  P  Drag  Coefficient  Lift  Coefficient 


d  Error 

Average 

Amplitud 

e 

Period 

(s) 

Averag 

e 

Amplitu 

de 

Period 

(s) 

10% 

2 

1.386 

0.195 

2.389 

0.001 

1.340 

4.768 

10% 

3 

1.370 

0.199 

2.369 

0.002 

1.327 

4.743 

10% 

4 

1.375 

0.202 

2.368 

0.005 

1.339 

4.743 

10% 

5 

1.375 

0.202 

2.372 

0.006 

1.340 

4.741 

5% 

3 

1.422 

0.216 

2.352 

-0.002 

1.402 

4.723 

5% 

4 

1.422 

0.217 

2.357 

0.004 

1.404 

4.720 

5% 

5 

1.420 

0.217 

2.359 

0.004 

1.405 

4.727 
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Figure  10.  Average  of  the  drag  coefficients  for  different  cases:  (a)  drag  coefficients; 


(b)  Lift  coefficients 


coefficients;  (b)  Lift  coefficients 
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(a)  (b) 

Figure  12.  Periods  of  the  drag  and  lift  coefficients  for  different  cases;  (a)  drag 

coefficients;  (b)  Lift  coefficients 

The  averages,  amplitudes,  and  periods  of  the  drag  and  lift  force  coefficients 
computed  from  the  Giid-II  and  Gird-III  are  plotted  in  Figures  10,  11,  and  12.  The 
results  obtained  from  Grid-I  are  also  plotted  in  these  figures  for  comparison 
purpose.  In  these  figures,  the  notations  of  the  six  cases  are  listed  in  Table  4. 
Obviously,  when  the  highest  polynomial  order  P  is  set  as  3,  the  accuracy  of  the 
results  obtained  from  Grid-II  and  Grid-Ill  can  improved  some.  However,  the  changes 
of  the  accuracy  of  the  drag  and  lift  coefficients  are  very  insignificant  for  the  P=4  and 
5.  Consequently,  P=4  for  the  present  case  is  usually  enough. 


Table  4.  Notation  of  the  six  cases 


Case 

Prescribed 

Error 

Grid 

Case 

Prescribed 

Error 

Grid 

1 

10% 

Grid-I 

4 

5% 

Grid-II 

2 

5% 

Grid-I 

5 

10% 

Grid-III 

3 

10% 

Grid-II 

6 

5% 

Grid-Ill 

The  maximum  errors  of  the  velocity  computed  between  50-second  and  60-second 
are  plotted  in  Figure  13.  Obviously,  the  errors  for  the  P=2  are  very  big  and  most  of 
them  are  higher  than  30%  for  Grid-I.  When  P  is  set  as  three,  the  accuracy  of  the 
velocity  improves  a  lot.  Unfortunately,  the  errors  do  not  change  much  when  P  is 
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increase  from  three  to  four  and  then  five.  There  are  two  possibilities.  One  is  that  the 
error  estimation  defined  in  this  paper  is  unreasonable.  The  other  is  that  higher 
order  polynomials,  4th  and  5th  for  example,  do  not  have  much  effect  on  the 
accuracy  of  velocity.  Sometimes,  the  errors  reduce  with  the  increase  of  the  highest 
polynomial  order  P.  The  errors  in  plot  13(b)  for  P=4  and  5,  for  example,  are  higher 
than  those  resulted  from  P=3.  One  possible  reason  is  that  the  orders  of  the 
interpolation  polynomials  are  set  as  2  for  all  elements. 

It  can  be  seen  from  plot  13(b)  that  the  errors  of  the  velocity  obtained  from  P=4  are 
much  higher  those  from  P=3.  However,  the  drag  and  lift  force  coefficients  for  both 
cases  are  very  close  that  are  indicated  in  Case  2  of  plots  10(a)  and  11.  The  same 
phenomenon  happened  for  Grid-II.  The  accuracy  of  the  velocity  for  P=4  and  5  is 
much  lower  than  that  for  P=3  while  the  aerodynamic  force  coefficients  are  very 
close.  If  the  error  estimation  defined  in  this  report  is  reasonable,  the  aerodynamic 
force  coefficients  are  insensitive  to  the  accuracy  of  the  velocity.  This  means  coarse 
grids  and  low  order  interpolation  polynomial  can  produce  similar  results  with 
respect  to  the  refined  grids  and  high  order  polynomial. 


(a)  (b) 

Figure  13.  Maximum  errors  between  50-second  and  60-second  for  Grid-I: 
(a)  Error  tolerance  10%;  Error  tolerance  5% 
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(a)  (b) 

Figure  14.  Maximum  errors  between  50-second  and  60-second  for  Grid-II: 
(a)  Error  tolerance  10%;  Error  tolerance  5% 


(a)  (b) 

Figure  15.  Maximum  errors  between  50-second  and  60-second  for  Grid-Ill: 

(a)  Error  tolerance  10%;  Error  tolerance  5% 

Together  with  the  numbers  of  equations  of  velocity  and  pressure  at  the  time  60- 
second,  the  run  times  for  the  three  grids  are  listed  in  Tables  5,  6,  and  7  for 
difference  Ps.  The  run  times  mentioned  here  are  not  the  exact  and  can  only  be  used 
qualitatively. 

Let  us  look  at  the  results  for  Grid-I  with  the  error  tolerance  is  10%  that  are  listed  in 
the  up  part  of  Table  5.  When  the  P  increases  from  2  to  3,  the  number  of  equations 
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for  velocity  and  pressure  increase  by  16.9%  and  292%  respectively.  Unfortunately, 
the  drag  force  coefficient  reduces  a  lot  which  leads  to  a  much  higher  error  that  P=2 
because  the  “exact”  coefficient  for  this  grid  is  looked  as  1.42.  After  that,  even 
though  the  coefficient  increase  when  P  increases  to  4  and  5,  they  are  still  lower 
than  that  for  P=2.  Hence,  for  this  case  the  p-adaptive  technique  is  absolutely 
unnecessary.  For  the  case  with  error  tolerance  5%,  the  result  seems  a  little  better. 
When  P  increases  from  2  to  3,  the  number  of  equations  of  velocity  and  pressure 
increase  by  28.3%  and  293%  respectively.  The  accuracy  of  the  corresponding 
aerodynamic  coefficients  also  increases  significantly.  Unfortunately,  even  though 
the  number  of  equations  increases  a  lot,  the  drag  and  lift  coefficients  change  very 
insignificant.  The  similar  phenomenon  happens  in  the  Grid-II  and  Grid-Ill. 


Table  5.  Number  of  equations  and  run  time  for  Grid-I 


Error 

Pol5momial 

Order 

Eqs.  Of  Velocity  At 

60s 

Eqs.  of  Pressure  At 

60s 

Run  Time 

(min) 

10% 

2 

8118 

2064 

209 

10% 

3 

9492 

8094 

421 

10% 

4 

9758 

9613 

459 

10% 

5 

10524 

10447 

670 

5% 

3 

10414 

8118 

554 

5% 

4 

10972 

10076 

661 

5% 

5 

14579 

14446 

1855 

Table  6.  Number  of  equations  and  run  time  for  Grid-II 


Error 

Polynomial 

Order 

Eqs.  Of  Velocity  At 

60s 

Eqs.  of  Pressure  At 

60s 

Run  Time 

(min) 

10% 

2 

12948 

3279 

365 

10% 

3 

14291 

12909 

630 

10% 

4 

14334 

14126 

638 

10% 

5 

14398 

14236 

646 

5% 

3 

15994 

12948 

899 

5% 

4 

16220 

15527 

881 

5% 

5 

17693 

17522 

1368 
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Table  7.  Number  of  equations  and  run  time  for  Grid-Ill 


Error 

Polynomial 

Order 

Eqs.  Of  Velocity  At 

60s 

Eqs.  of  Pressure  At 

60s 

Run  Time 

(min) 

2 

29610 

7470 

1125 

3 

30735 

29511 

1346 

4 

30742 

30511 

1405 

10% 

5 

30781 

30556 

1338 

5% 

3 

33698 

29610 

1889 

5% 

4 

33770 

33356 

2052 

5% 

5 

34084 

33801 

1948 

Effects  of  the  Error  Tolerance 

The  program  is  run  again  for  the  Gird-I  when  the  error  tolerances  are  set  as  2%  and 
1%  respectively.  The  maximum  errors  computed  between  50-second  and  60-second 
with  P=5  are  plotted  in  Figure  16.  The  percentages  of  the  elements  with  2nd,  3rd, 
4th,  and  5th  orders  interpolation  polynomial  are  listed  in  Table  8.  They  are  resulted 
from  P=5  and  error  tolerances  are  10%,  5%,  2%,  and  1%  respectively.  The 
corresponding  averages,  amplitudes,  and  periods  for  the  error  tolerances  2%  and 
1%  are  1.421,  0.217,  2.351  and  1.420,  0.215,  2.359  respectively. 

The  averaged  maximum  errors  between  50-second  and  60-second  are,  respectively, 
0.1102,  0.0885,  0.0702,  and  0.0760  for  t)=10%,  5%,  2%,  and  1%.  The  accuracy  of 

the  results  improves  a  lot  when  the  prescribed  error  changes  from  10%  to  5%. 
There  is  no  much  change  when  the  error  tolerance  reduces  from  5%  to  2%  and  then 
to  1%  even  though  the  percentage  of  the  elements  with  5th  order  polynomial 
increases  from  10.93%  to  20.40%  and  then  27.07%.  Again,  it  seems  higher  order 
polynomial  does  not  have  much  effect  on  the  accuracy  of  the  force  coefficients. 

As  for  the  aerodynamic  forces,  there  is  no  much  difference  for  the  Grid-II  amd  Grid- 
Ill  when  the  error  tolerance  reduces.  One  reason  is  that  the  results  obtained  from 
Grid-II  and  Grid-Ill  are  already  very  close  to  the  exact.  Again,  the  higher  order 
polynomial  seems  insignificant  to  the  accuracy  of  the  drag  and  lift  force  coefficients 
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for  the  present  problem.  Consequently,  although  the  low  error  tolerance  can 
increase  the  number  of  equations  of  velocity  and  pressure  and  the  accuracy  of  the 
velocity  somewhat,  it  can  not  increase  the  accuracy  of  the  aerodynamic  forces 
significantly. 


Figure  16.  Maximum  errors  between  50-second  and  60-second  for  different  error 

tolerances  with  P=5 


Table  8.  Percentages  of  the  elements  with  2nd,  3rd,  4th,  and  5th  orders  of 

polynomial 


Error  Tolerance 

Order  2 

Order  3 

Order  4 

Order  5 

10% 

83.26 

11.78 

3.21 

1.75 

5% 

71.88 

11.33 

5.86 

10.93 

2% 

62.01 

12.58 

5.01 

20.40 

1% 

59.00 

8.07 

5.86 

27.07 

In  the  above  discussion,  the  error  distribution  is  evaluated  every  one  second.  To 
make  sure  that  the  refinement  has  been  done  completely  within  the  60  seconds,  0.5 
second  and  0.2  second  are  used  for  the  time  interval.  The  computed  drag  and  lift 
coefficients  are  plotted  in  Figure  17.  Except  some  time  shift,  the  amplitudes  and 
averages  of  the  drag  and  lift  coefficients  are  very  close.  This  means  that  results 
provided  above  are  reasonable. 
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(a)  (b) 

Figure  17.  Drag  and  lift  coefficients  for  different  time  interval  of  refinement; 
(a)  Drag  coefficients;  (b)  List  Coefficients 


Effects  of  the  Grid 

It  can  be  seen  from  Figures  13(b),  14(b),  and  15(b)  that  the  cases  of  P=3,  P=3,  and 
P=3  or  4  or  5  have  the  highest  accuracy  of  velocity  for  Grid-I,  Grid-II,  and  Grid-Ill 
respectively.  Hence,  the  aerodynamic  forces  of  these  three  cases,  which  are  listed  in 
Table  9,  are  considered  as  the  exact  for  each  grid.  For  the  Grid-I  with  P=3,  the 
number  of  equations  of  velocity  and  pressure  are  10414  and  8118.  They  are  15994, 
12948  and  33698,  29610  for  Grid-II  and  Grid-Ill  with  P=3  respectively.  The 
differences  of  the  aerodynamic  forces  among  the  three  grids  are  very  minor  while 
the  differences  of  the  number  of  equations  are  very  significant.  Consequently,  it  is 
unnecessary  to  use  highly  refined  grid  to  evaluate  the  drag  and  lift  forces. 


Table  9.  The  exact  aerodynamic  forces  for  the  three  grids 


Grid 

Average  Cd 

Amplitude 

Cd 

Period  (s)  Cd 

Amplitude 

Cl 

Period  (s)  Cl 

I 

1.422 

0.216 

2.352 

1.402 

4.723 

II 

1.428 

0.223 

2.320 

1.414 

4.646 

III 

1.444 

0.226 

2.309 

1.439 

4.620 
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S.Conclusions 

The  p-version  adaptive  finite  element  method  has  been  implemented  into  the 
computational  flow  around  a  circular  cylinder  to  compute  aerodynamic  forces.  The 
second  through  fifth  orders  of  polynomials  are  considered  for  the  velocity  in  the 
adaptive  method.  One  order  less  of  pol5momial  is  used  for  pressure.  Velocity  is 
selected  as  the  error  estimator.  A  flow  around  a  circular  cylinder  with  Reynolds 
number  of  1000  was  simulated  using  this  technique.  The  effects  of  the  highest 
order  of  polynomials,  error  tolerance,  and  size  of  the  element  on  the  accuracy  of  the 
drag  and  lift  coefficients  are  surveyed  using  this  flow  simulation. 

When  the  grid  is  coeirse,  the  p-adaptive  technique  is  very  efficient.  It  can  improve 
accuracy  significantly  while  the  computed  time  does  not  increase  much.  For  the 
refined  grid,  the  results  of  the  present  simulation  show  that  the  higher  orders  of 
interpolation  polynomials,  P=4  and  5  for  example,  do  not  have  much  effect  on  the 
accuracy  of  the  velocity  as  well  as  the  drag  and  lift  coefficients.  The  error  of  the 
velocity  does  not  have  much  effect  on  the  accuracy  of  the  drag  and  lift  force 
coefficients.  If  the  error  estimation  defined  in  this  report  is  reasonable,  the 
aerodynamic  force  coefficients  are  insensitive  to  the  accuracy  of  the  velocity.  The 
accuracy  of  the  results  improves  a  lot  when  the  prescribed  error  changes  from  10% 
to  5%.  There  is  no  much  change  when  the  error  tolerance  reduces  from  5%  to  2% 
and  then  to  1%  even  though  the  percentage  of  the  elements  with  5th  order 
polynomial  increases.  The  differences  of  the  aerodynamic  forces  among  the  three 
refined  grids  are  very  minor  while  the  differences  of  the  number  of  equations  are 
very  significant. 

The  code  will  be  applied  to  the  practical  aerospace  problems,  especially  the 
AFOSR/DEPSCoR  funded  work  on  fluid-structure  interaction  of  aerospace 
structures. 
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1.  Introduction 

The  development  of  the  advancing  front  mesh  generation  method  is  described  in 
a  series  of  papers  by  Peraire  et  al  [1-4].  Important  contributions  to  the  development 
of  the  method  have  also  been  given  by  Lohner  et  al  [5-7].  Although,  two  dimensional 
mesh  generation  methods,  similar  in  concept  to  the  advancing  front  method,  have 
been  used  since  at  least  the  early  seventies,  the  generalizations  which  make  the 
advancing  front  method  qualify  as  an  automatic  three  dimensional  mesh  generation 
method  were  not  devised  until  the  mid-eighties  by  the  previously  referenced  works 
by  Peraire  et  al  [3,4]. 

The  characteristic  feature  of  the  advancing  front  method  is  that  [7]  elements  and 
nodes  axe  created  simultaneously,  with  all  elements  residing  behind  a  front  which 
sweeps  over  the  domain,  starting  from  external  as  well  as  internal  boundaries.  The 
front  consists  of  those  mesh  line  segments  that  separates  the  discreth^d  and 
undiscretized  parts  of  the  domain  from  each  other.  Hence,  the  initial  front  consists 
of  the  mesh  line  segments  that  make  up  the  boundary  discretization.  A  line 
segment  is  selected  from  the  front  and  a  new  node  is  constructed  so  that  a  new 
triangular  element  is  formed.  As  this  procedure  is  repeated,  the  front  advances  over 
the  domain  and  changes  continuously  in  that  old  line  segments  are  deleted  and 
new  added  as  triangles  are  created.  The  mesh  is  complete  when  the  front  is  empty. 


2.  Generation  of  Initial  Mesh 

2.1.  Background  Grid 

The  problem  of  generating  a  mesh  over  a  two  dimensional  region  of  arbitrary 
shape  is  considerably  simplified  if  unstructured  triangular  meshes  are  employed. 
For  the  method  to  be  described  here,  this  process  is  started  by  constructing  by 
hand  a  coarse  background  grid  of  3-node  triangular  elements  which  completely 
covers  the  solution  domain  of  interest  [2].  The  background  grid  is  used  to  provide  a 
piecewise  linear  spatial  distribution  for  these  parameters  over  the  grid  to  be 
generated.  Thus,  at  each  node  on  the  background  grid,  the  stretching:  node  spacing 
6  ,  value  of  stretching  parameter  s  and  the  direction  of  stretching  a  must  be 
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specified.  During  the  generation  process  the  local  values  of  these  quantities  will  be 
obtained  by  linear  interpolation,  over  the  triangles  of  the  background  grid,  between 
the  specified  nodal  values.  If  5  is  required  to  be  uniform  initially  and  no  stretching 
is  to  be  specified,  then  the  background  grid  need  only  consist  of  a  single  element 
which  covers  the  solution  domain. 

2.2.  Generation  of  Boundary  Nodes 

Define  the  boundaries  of  the  domain  to  be  gridded.  This  is  typically 
accomplished  by  spines  in  2-D  and  surface  patches  in  3-D.  The  boundary  of  the 
solution  domain  is  represented  by  the  union  of  closed  loops  of  curved  segments  and 
boundary  nodes  are  placed  at  the  points  of  intersection  of  these  segments.  For 
simply  connected  regions  there  is  only  one  closed  loop,  whereas  for  multi-connected 
regions  there  will  be  as  many  internal  loops  as  the  number  of  openings  inside  the 
domain.  The  segments  of  exterior  boundary  are  defined  in  an  anti-clockwise 
manner  while  the  segments  of  the  interior  boundaries  are  specified  in  a  clockwise 
fashion.  Tliis  means  that,  as  the  boundary  curve  is  traversed,  the  region  to  be 
triangulated  always  lies  to  the  left.  Before  beginning  the  process  of  generating 
triangles  within  the  region  of  interest,  the  positioning  of  additional  nodes  on  the 
boundaries  of  the  region  has  to  be  performed.  Each  boundary  segment  is 
considered  in  turn  and  nodal  points  are  generated  on  the  boundary  segments,  with 
the  spacing  of  the  points  being  determined  by  interpolated  values  of  5  ,  J  and  a 
[2].  This  yields  the  initial  front. 

2.3.  Triangle  Generation 

At  the  start  of  the  process  the  front  consists  of  the  sequence  of  straight  line 
segments  which  connect  consecutive  boundary  nodes.  During  the  generation 
process,  any  straight  line  segment  which  is  available  to  form  an  element  side  is 
termed  active,  whereas  any  segment  which  is  no  longer  active  is  removed  from  the 
front.  Thus,  while  the  domadn  boundary  will  always  remain  the  same,  the 
generation  front  will  change  continuously  and  has  to  be  updated  whenever  a  new 
element  is  formed.  The  following  steps  are  involved  in  the  process  of  generation  a 
new  triangle  in  the  mesh. 
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2.3.1  Prepare  the  creation  the  new  element 

2.3. 1.1  Select  the  base  edge  on  which  the  element  is  to  be  constructed.  If  large 
variations  in  6  are  present  in  the  background  grid  then  it  is  advantageous 
to  look  for  the  smallest  active  side,  but  if  8  is  constant  or  varying  slowly 
the  last  active  side  in  the  front  is  used. 

2.3. 1.2  Determine  the  element  size,  i.e.  the  length  of  the  two  new  edges  of  the  so- 
called  ideal  element.  Suppose  the  chosen  side  joins  nodes  A  and  B. 
Determine  the  local  mesh  parameters  5^ ,  and  at  the  mid-point  of 
AB  by  interpolating  over  the  background  grid. 

^  ‘  ^MFG  ‘  ^MGE  '  ^MEF  JJJ 

^EFG 

where  Sj^p^ ,  8,^^^ ,  S^^^p  and  Spp^  are  the  areas  of  the  triangles  MFG,MGE, 
MEF  and  EFG  and  they  are  defined  in  Figure  1.  Make  a  local  rotation  of 
coordinates  so  that  lies  along  the  x,  axis  and  scale  the  coordinate  by 

a  factor  s^f . 

2.3. 1.3  Define  the  spatial  neighborhood,  i.e.  the  region  close  to  the  base  edge  where 
the  configuration  of  the  front  must  be  known  to  make  the  validation  of  the 
new  element  possible. 


Figure  1.  Area  distribution 


2.3.2  Create  the  new  element 

2.3.2. 1  Determine  the  position  of  the  ideal  node  so  that  the  resulting  element  is  as 
equilateral  as  possible  while  respecting  the  element  size  as  defined  by 
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background  grid.  In  the  new  coordinate  system,  a  triangle  which  is  as 
regular  as  possible  will  be  generated.  Determine  5^  according  to 

'0.55AB  5^<0.55AB 

0.55  AB  <5^<  2AB  (2) 

2AB  2AB  < 

The  inequalities  used  here  are  necessary  to  obtain  ecometrical  compatibility 
and  to  ensure  that  the  elements  with  essive  distoration  are  not  generated. 
Different  inequalities  can  be  devised  out  values  shown  have  worked  well  in 
practice.  Note  that  not  only  the  element  size,  but  also  the  configuration  of 
the  front  is  .considered  when  determining  the  position  of  the  ideal  node. 
Also  note  that  the  configuration  of  the  front  might  prevent  creation  of  the 
ideal  node.  If  no  node  created  go  to  2.3. 2.3. 

2. 3. 2. 2  Check  the  validity  and  the  suitability  of  the  ideal  element.  The  validity 
check  consists  in  ensuring  that  none  of  the  newly  created  edges  intersects 
with  the  front  in  the  neighborhood.  The  suitability  check  consists  in 
ensuring  that  the  node  and  the  edges  of  the  element  are  sufficient  far  from 
nodes  and  edges  in  the  front.  If  the  ideal  element  satisfies  the  validity  and 
suitability  requirements  go  to  2.3.3. 

2. 3. 2. 3  Check  if  any  existing  active  node  in  the  neighborhood  can  be  used  to  form 
the  new  element.  A  node  should  be  positioned  sufficiently  close  to  the  base 
edge  to  be  considered  a  candidate  node.  Usually,  the  candidate  nodes  are 
defined  as  those  lie  within  the  circle  with  centre  at  C  and  radius  nAB. 
(There  is  no  unique  choice  for  the  value  of  n  which  should  be  adpoted,  but 
the  value  2  has  been  used  for  the  program.)  These  nodes  are  odered 
according  to  their  distance  from  C  with  the  first  node  in  the  list  being  the 
closest  to  C.  The  validity  and  the  suitability  of  the  elements  formed  with  the 
candidate  nodes  are  checked  in  accordance  with  step  2. 3. 2. 2.  If  several 
elements  pass  both  these  checks,  select  the  best  element  as  determined  by 
the  suitability  check  and  then  go  to  2.3.3. 

2. 3. 2.4  Create  a  set  of  try  nodes  on  positions  which  are  predetermined  relative  the 
base  edge.  The  vadidity  and  the  suitability  of  the  elements  formed  with  the 
try  nodes  are  checked  in  accordance  with  step  2. 3. 2.2.  If  several  elements 
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pass  both  these  checks,  select  the  best  element  as  determined  by  the 
suitability  check  and  then  go  to  2.3.3. 

2. 3. 2. 5  Check  the  validity  and  the  suitability  of  the  elements  formed  with  the 
existing  nodes  in  the  neighborhood  as  well  as  with  the  try  nodes  in  step 
2.3. 2.4.  The  The  validity  and  the  suitability  of  the  elements  formed  with  the 
try  nodes  nodes  are  checked  in  accordance  with  step  2. 3. 2. 2.  Among  those 
elements  that  pass  the  validity  check,  select  the  best  element  as 
determined  by  the  suitability  check  and  then  go  to  2.3.3. 

2.3.3  Update  the  front  data 

First,  delete  edges  and  nodes  that  no  longer  are  members  of  the  front.  Second, 
add  a  new  node  and  new  edges  to  the  front  data  structures.  The,  if  the  front  is  not 
empty  go  to  2.3.2. 1  otherwise  finish. 

2.4  Implementations 
2.4.1  Normalized  Space 

The  construction  of  the  ideal  triangle  is  done  in  a  transformed,  so-called 
normalized  space  as  proposed  by  Peraire  et  al  [2].  This  is  not  absolutely  necessary 
but  convenient  in  case  of  anisotropic  mesh  control,  i.e.  when  the  size  of  the  triangle 
should  be  different  in  mutually  orthogonal  directions.  Anisotropic  mesh  control  is 
most  easily  specified  by  a  node  spacing  6  ,  a  value  of  stretching  parameter  s  and  a 
direction  of  stretching  a  .  Note  that  we  assume  a  to  be  a  unit  vector.  The  variation 
over  the  domain  of  the  mesh  characteristics  is  given  by  a  mesh  size  function,  for 
instance,  a  background  mesh. 

The  first  step  of  the  coordinate  transformation  from  model  space  to  normalized 
space  is  a  translation  of  the  origin  of  the  model  space  coordinate  system  to  the 
midpoint  of  the  base  edge,  .  In  the  second  step,  the  coordinate  system  is  rotated 
so  that  the  axis  is  aligned  with  the  stretch  direction  a  .  The  third  step  is  to  scale 
with  l/s  along  the  a  direction  by  which  a  normailized  space  with  isotropic  mesh 

control  is  obtained.  Hence,  the  ideal  triangle  is  equilateral  and  the  element  size  is  6 
in  the  normalized  space.  The  coordinate  transformation  from  model  space  to 
normalized  space  is  most  easily  expressed  as 
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N  =  SRT  = 
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where  the  transformation  matrix  N  is  composed  of  the  translation  matrix  T ,  the 
rotation  matrix  R  and  the  scaling  matrix  S .  The  mesh  parameters  s ,  5  and  a 
defining  the  transformation  are,  for  instance,  evaluated  at  the  midpoint  of  the  base 
edge.  0  is  defined  the  angle  from  to  x^  and 
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Note  that  the  transformation  matrix  defined  by  equation  ()  assumes  homogeneous 
coordinates  which  iin  this  context  simply  means  that  a  two  dimensional  point 

jc  =  [x,  is  expressed  as  the  triple  =  [x,  X2  l]^.  Furthermore,  note  that  the 

row  vectors  of  R  rotate  into  the  new  coordinate  axes.  The  inverse  transformation 
from  the  normalized  space  to  model  space  is 


M  =  rt-'5"'  = 


'1 

0 

•^m.l 

COS  6 

-sin^ 

O' 

s 

0 

O' 

0 

1 

^m,2 

sin^ 

cos^ 

0 

0 

1 

0 

0 

0 

1  _ 

0 

0 

1 

_0 

0 

1 

(5) 


2.4.2  Validity  Check 

Candidate  nodes  and  trian^es  must  satisfy  certain  so-called  validity  criteria  to 
be  accepted  as  member  in  the  mesh.  The  criterion  whether  to  accept  a  candidate 
node  is  that  it  should  be  positioned  sufficiently  far  from  existing  nodes  on  the  front. 
Here  sufficiently  far  is  0.67J^j  which  is  the  element  size.  Hence,  the  node  and 
triangle  validity  criteria  are  as  follows: 

•  A  new  node  is  valid  and  can  be  accepted  as  a  member  of  the  mesh  if  the 
distance  from  it  to  the  closest  node  is  larger  than  0.615^^ . 

•  A  triangle  is  vald  and  can  be  accepted  as  a  member  of  the  mesh  if  the  two  edges 
do  not  intersect  the  advancing  front. 
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The  node  validity  criterion  is  simply  that  a  new  candidate  node  should  be  at 
least  the  distance  0.67^^^  away  from  existing  nodes  on  the  front.  To  check  the 

validity  of  a  node  is  straightforward  since  it  only  consists  of  distance  calculations. 
For  efficiency  reasons  square  root  evaluations  should  be  avoided.  Another 
optimization  is  to  check  the  distance  in  each  coordinate  direction  separately  in  a 
filtering  step  by  which  nodes  out  of  bounds  are  quickly  discarded.  Note  that  the 
neighborhood  is  large  enough  to  contain  all  nodes  of  interest  for  the  node  validity 
check. 

The  element  validity  criterion  is  what  a  new  element  is  valid  if  it  does  not 
intersect  the  front.  To  ensure  the  validity  of  an  element  is  thus  equivalent  with  to 
check  if  any  front  edge  intersects  any  of  the  edges  of  the  candidate  element.  Here 
we  assume  that  the  elements  are  straight  sided  triangles,  and  the  validity  check 
then  reduce  to  determine  whether  bounded  line  segment  intersects.  Even  though 
this  is  rather  straightforward,  the  great  number  of  intersections  calculations  that 
needs  to  be  done  for  the  construction  of  a  mesh  motivates  a  carefully  designed 
algorithm.  The  core  of  the  algorithm  we  propose  is  based  on  a  line  clipping 
algorithm  described  by  Blinn  [3].  Before  we  proceed  with  description  of  the 
algorithm  for  determining  whether  a  front  edge  intersects  a  triangle  in  two 
dimensional  space,  we  introduce  some  useful  notation.  Each  triangle  edge,  E-^ , 
1=1, 2, 3  define  a  line  that  divides  the  Euclidian  plane  into  two  half  spaces  which  we 
call  the  positive  half  space,  and  the  negative  half  space,  Ef_.  The  triangle  is 

oriented  so  that  its  interior  is  in  ,  and  the  normal  vector  ,  tj  ,  to  £',■  is  directed 


towards  the  interior  of  the  triangle.  The  implicit  equation  of  the  line  E.  is 


(6) 


where  is  a  point  on  the  line,  x  is  an  arbitray  point,  and  b^‘  is  the  so  called 
boundary  coordinate  of  the  point  x  with  respect  to  the  line  Ei .  The  following 
relations  holds  for  the  boundary  coordinates 

f>0 


=  0  (x^E,) 

<0  (xeE,.) 


(7) 
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Note  that  lb: 


equals  the  perpendicular  distance  between  the  point  x  and  the  line 


Ej  if  n  is  a  unit  vector.  Also  note  that  we  use  the  term  line  and  edge  to  distinguish 


between  the  unbounded  £ind  the  bounded  line,  respectively. 

A  three  layered  filtering  algorithm  for  checking  whether  a  candidate  triangle,  t , 
intersects  a  front  edge,  e ,  is  defined  as  follows: 

1.  Check  if  both  endpoints  of  e  are  contained  in  the  one  and  same  negative  half 
space  as  defined  by  any  of  the  triangle  edges  .  If  that  is  no  intersection 
between  e  and  t,  otherwise  go  to  2. 

2.  Check  if  one  or  both  endpoints  of  e  are  in  the  interior  of  i .  If  that  is  the  case, 
there  is  an  intersection  between  e  and  t ,  otherwise  go  to  3. 

3.  Check  for  intersection  between  e  and  E^,  where  1=1,2.  If  there  is  no 
intersection,  it  can  be  finally  concluded  that  e  and  t  does  not  intersect. 

Step  1  is  based  on  the  observation  that  intersection  cein  be  ruled  out  if  both 
endpoints  of  the  front  edge  are  in  the  negative  half-space  of  at  least  one  of  the 
triangle  edges.  Note  that  the  boundary  coordinates  for  the  front  edge  are  saved  as 
they  are  calculated  since  they  also  are  used  in  steps  2  and  3.  Thus,  the  front  edge 
and  the  triangle  do  not  intersect  if  the  boundary  coordinates  for  both  endpoints  of 
the  edge  are  contained  in  the  triangle.  A  point  is  contained  in  the  triangle  if  and 
only  if  the  point  is  in  the  positive  half-space  for  each  one  of  the  triangle  edges.  This 
corresponds  to  that  aU  three  boundary  coordinates  for  a  front  edge  (one  for  each  E^ ) 


are  positive.  Step  3  is  reached  only  if  the  front  edge  intersects  at  least  one  of  the 
straight  lines  defined  by  the  trian^e  edges  E-  (not  necessarily  within  the 

boundaries  of  the  triangle  edge).  The  check  in  step  3  is  based  on  that  there  is  an 
intersection  between  a  line  segment  in  parametric  form 

a:  =  Xo +5(x, -Xq)  (8) 


and  a  triangle  edge  E-  if  and  only  if  the  boundary  coordinates  6^'  and  6^'  have 


opposite  signs,  i.e.  the  points  x,,  and  x,  are  on  different  sides  of  the  line  E^ .  The 
point  of  intersection  between  the  line  segment  and  the  line  is  then  obtained  from 
equation  (6)  and  (7),  viz.. 
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An  intersection  s  is  calculated  for  each  line  for  which  the  boundaiy  coordinates 

b^‘  and  have  opposite  signs.  What  remains,  is  to  determine  whether  any  of  the 

intersections  is  within  the  boundaries  of  a  triangle  edge.  This  is  accomplished  by 
associating  the  parameters  and  with  the  points  and  x^ ,  respectively.  The 

parameters  and  are  initialized  to  zero  and  one,  respectively.  For  each 

intersection  calculated  from  equation  (9),  the  parameter  whose  associated  boundary 
coordinate  is  negative  is  updated,  viz., 

k=max(/o,5)  (6j<0) 

[i,  =min(/„5)  (6^^' <o) 

If  /q  ever  becomes  greater  than  (/q  >^i)  the  line  segment  and  the  triangle  does  not 
intersect.  Consequently,  the  Line  segment  crosses  the  triangle  if  io  still  less  than 
<^,)  when  e  has  been  checked  against  aU  three  .  Finally,  note  that  the  second 

and  third  step  can  be  done  simultaneously  since  the  computations  done  in  step  two 
also  must  be  done  in  step  three  (check  for  boundary  coordinate  with  opposite 
signs).  This  is  usually  efficient  since  the  inclusion  of  front  nodes  (step  two)  in  a 
triangle  is  not  so  common. 

2.4.3  Suitability  Criteria 

The  satisfaction  of  the  validity  criteria  ensures  the  creation  of  a  topologically 
compatible  mesh.  However,  the  validity  criterion  does  not  ensure  the  creation  of 
well  shaped  elements.  To  avoid  the  construction  of  ill  shaped  elements,  it  is 
necessary  to  control  the  distance  between  nodes  and  edges  as  the  front  advances 
over  the  domain.  Based  on  the  considerations  above  we  try  to  create  well  behaved 
meshes  without  having  to  resort  to  global  smoothing.  To  achieve  this  goal  it  is 
necessary  to  control  the  distance  between  nodes  and  edges  before  accepting  any 
new  node  or  edge.  In  fact,  the  suitability  criteria  are  expressed  as  an  angle  criterion 
which  is  as  follows: 
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•  A  new  node  is  accepted  as  a  member  of  the  mesh  if  the  smallest  angle  formed  by 
the  node  and  any  front  edge  is  larger  than  >  30° . 

•  A  new  edge  is  accepted  as  a  member  of  the  mesh  if  the  smallest  angle  formed  by 
the  edge  and  any  front  is  >  30° . 


2.4.4  Control  Line  Discretization 

Curve  discretization  is  the  least  covered  topic  of  mesh  generation.  Perhaps 
because  the  subject  is  viewed  as  trivial.  Although,  curve  discretization  algorithms 
have  been  proposed  in  References  [3,6,9,10]  and  discussed  by  Frykestig  in  detail 


2.4.4. 1  Straight  Linear  Interpolation 

Straight  line  interpolation  is  the  simplest  case  for  the  discretization  of  boundary 
lines.  It  is  used  for  uniformed  element  size.  The  number  of  nodes  or  elements  that 
should  be  created  on  the  line  is  usually  specified  in  this  case.  The  nodal  locations 
are  determined  through  interpolation.  Usually  linear  interpolation  which  results  in 
an  equal  subdivision  of  the  line  is  used. 

When  the  number  of  elements,  ,  is  defined,  the  element  length  is 


I 


N, 


(11) 


where  is  the  total  length  of  the  boundary  line.  If  the  length  of  elements,  /  ,  is 
specified,  the  number  of  elements  is 

(12) 


N,= 


f^O.5 


Int 


Then  equation  (11)  is  used  to  calculate  the  actual  element  length.  The  locations  of 
node  «,•  is 


(13) 


2.4.4.2  Arithmetic  series  expansion 

The  use  of  arithmetic  [11]  and  geometric  [10]  series  expansions  is  quite  common 
for  line  discretization.  In  these  approaches,  it  is  required  to  specify  the  element  size 
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at  one  or  both  endpoints  of  the  line,  and,  possibly,  the  number  of  elements  to  be 
created  on  the  line.  These  methods  are  usually  used  in  the  case  that  the  element 
size  changes  slowly. 

Consider  a  line  with  prescribed  element  lengths  /,  and  I2  at  the  two  endpoints. 
They  can  be  obtained  from  the  background  definition.  An  arithmetic  series  line 
discretization  is  defined  as 

l,^l,.,+d  =  ll+{i-\)d  (14) 


,  NXi:+ll)  Af.f2/;  +  (W.-lV 


where  is  a  constant  for  the  line  that  determines  the  mesh  grading.  The  number  of 
elements  on  the  line  is  solved  from  equation  (15)  as 


K-  ^+0.5 

L'l  ^  h  J/„, 


The  lengths  of  the  elements  are  then  obtained  as 

/  +  /=i,2,.-.,A^,  (17) 

We  note  that  calculated  from  equation  (17),  is  likely  to  differ  from  the  specified 
element  size  /* .  The  residual  is  r  =  /*  -  .  It  will  satisfy  1^^  =  by  distributing  r  to 

the  interior  elements,  i.e. 


2.4.4.3  Meshing  size  function  approach 

In  this  method,  the  spacing  function  6(j)  is  defined  over  the  line  segment.  It 
defines  the  element  size,  or  the  distance  between  two  neighboring  nodes  at  the 
point  given  by  r(5)  equivalently.  Since  6  (s)  can  be  thought  on  as  having  the  unit 

lenght Inode,  the  reciprocal  \l8{s)  can  be  interpreted  as  a  node  density.  The  number  of 
element  which  need  to  be  created  along  the  boundary  line  is  calculated  by  direct 
integration  [3,9] 
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(19) 


Taking  equal  to  the  nearest  integer  to  A^.  The  nodal  positions,  s,. ,  are  calculated 
from 


(20) 


Suppose  the  element  sizes  at  both  endpoints  are  and  51  respectively.  The 
element  size  is  then  assumed  to  vary  linearly  between  the  two  endpoints.  The 
spacing  function  5  (5)  is 


Substituting  equation  (21)  into  equations  (19)  and  (20),  one  has 


Si  —  *^0  + 


O2 


(21) 


(22) 

(23) 


2.4.5  Locating  Query  Point  M 

In  the  h-version  adaptive  finite  element  method,  the  finite  element  gird  is 
generated  adaptively  according  to  the  error  in  the  computation.  Usually,  the 
adaptive  mesh  is  generated  by  considering  the  current  computational  gird  as  a 
background  gird  [2].  Because  the  number  of  elements  in  the  background  mesh  is 
large,  some  consideration  must  be  given  to  the  search  problem  of  finding  the 
element  in  the  background  mesh  in  which  the  query  point  M  is  located. 

2.4.5. 1  Searching  Algorithms 
Algorithm  I 

This  algorithm  requires,  for  each  element  e  of  the  background  grid,  the 
knowledge  of  the  three  surrounding  elements  which  have  sides  in  common  with 
element  e .  Given  the  coordinates  of  M  and  a  starting  element  of  the  background 
grid,  the  three  area-coordinates  of  M  are  determined.  If  each  area  coordinate  lies 
between  zero  and  one,  then  the  element  contains  the  point  M .  If  not,  the  node  for 
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which  the  area  coordinate  is  a  minimum  is  found  and  this  indicate  the  nest  element 
to  be  checked.  As  shown  in  Figure  2,  if  the  area  coordinates  Z, ,  and  of 

element  e  are  evaluated  at  M  and  i,  <  Z,2  <  Z-j ,  the  next  element  to  be  checked  is 

element  e, .  In  this  manner,  the  necessary  of  searching  over  all  the  elements  in  the 
background  grid  is  avoided. 


M 


Figure  2.  Searching  algorithm  to  locate  point  M  on  the  background  grid 


Algorithm  II 

Another  possibility  is  to  store  the  nodes  of  the  background  mesh  in  an  actree  as 
proposed  by  Lohner  [5].  Every  time  the  mesh  size  is  queried  at  some  point,  the  node 
of  the  background  mesh  closest  to  the  query  point  is  found  by  searching  the  octree. 
Once  the  closest  point  has  been  found,  it  is  straightforward  to  find  the  element 
containing  the  query  point  assuming  that  the  adjacency  relationships  between 
nodes  and  elements  in  the  background  mesh  are  available.  The  details  of  the  octree 
may  be  found  in  References  5,8. 


2.4.5.2  Locating  the  Point  M 
Algorithm  I  -  Peraire:  [2] 

If  the  three  area-coordinates  of  M  with  respect  to  a  triangle  element  in  the 
background  grid  lie  between  zero  and  one,  the  element  contains  the  point.  This 
condition  is  equivalent  to  that  the  smallest  area  coordinate  is  non-negative.  It  can 
be  proven  simply.  Suppose  coordinate  L-  is  the  smallest  among  the  three 


43 


Paer  II:  H-Adaptive  Advancing  Grid  Generation 


Selvam  and  Qu 


coordinates.  Because  Z,  >  0  and  L-  +  Lj  +  L,^  =  1 ,  we  have  Lj  +  <  \ .  With 

considering  Lj  >  L,.  >  0  and  >  L,  >  0 ,  the  conclusion  may  be  obtained. 

Algorithm  II  -  Frykestig  [8:  151-152] 

Each  triangle  edge,  ,  1=1, 2, 3  define  a  line  that  divides  the  Euclidian  plane  into 
two  half  spaces  which  we  call  the  positive  half  space,  and  the  negative  half 
space,  E._.  The  triangle  is  oriented  so  that  its  interior  is  in  E;^,  and  the  normal 
vector,  n  ,  to  E^  is  directed  towards  the  interior  of  the  triangle.  The  implicit  equation 
of  the  line  E.  is 

b^‘=n^{x-x,)  (24) 

where  Xg_  is  a  point  on  the  line,  x  is  an  arbitrary  point,  and  b^‘  is  the  so  called 
boundary  coordinate  of  the  point  x  with  respect  to  the  line  £) .  The  following 
relation  holds  for  the  boundary  coordinates 
>0  (xeEiJ 

1=0  (x  g  El)  (25) 

<0  {xe  Ei_) 

Note  that  |  equals  the  perpendicular  distance  between  the  point  x  and  the  line 

El  if  n  is  a  unit  vector.  Also  note  that  we  use  the  term  line  and  edge  to  distinguish 

between  the  unbounded  and  the  bounded  line,  respectively. 

For  the  line  segment  AB,  one  of  its  normal  vectors  n  through  the  original  may 
be  expressed  as 

n :  (0,0)  ^  (y^  -ys^x^-x^)  (26) 

Therefore,  equation  (24)  becomes 

=iyA-  yaXx  -  xe)+{xs  -x^Xy-  ys)  (27) 


3.  Adaptive  Remeshm< 


The  procedure  outlined  above  enable  an  initial  approximation  to  the  steady  state 
solution  to  be  obtained  for  a  given  problem.  The  solution  quality  can  be  improved  by 
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adaptively  refining  the  mesh.  This  mesh  adaptation  is  achieved  by  using  the 
computed  solution  to  determine  "optimum"  nodal  values  for  b  ,  s  ,  and  a .  The  mesh 
is  then  regenerated  with  the  initial  computational  mesh  now  acting  as  a 
background  grid. 


4.  Numerical  Examples 

Example  I 

As  shown  in  Figure  3,  the  size  of  the  rectangle  ABCD  is  Ix2xlx2.  Four 
triangular  elements  are  used  for  the  background  grid.  At  first,  the  uniformed  grid 
size  is  used.  Therefore,  the  input  data  are  given  by 


A 

4 


5 

B 


D 

3 


2 

C 


Figure  3.  Rectangle  to  be  meshed 


INODE 

NISEG 

NBACK 

NPOIN 

4 

4 

4 

5 

Coordinates  of  the  boundary  nodes: 


X 

Y 

0 

1 

0 

0 

2 

0 

2 

1 

Line  segments  defining  the  boundary: 
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2 

3 

3 

4 

4 

1 

Definition  of  the  grid  size  through  background  element: 


X 

Y 

Size 

1.5 

0.5 

0.1 

2 

0 

0.1 

2 

1 

0.1 

0 

1 

0.1 

0 

0 

0.1 

Connectivity  of  the  background  grid: 


1 

2 

3 

1 

3 

4 

1 

4 

5 

1 

5 

2 

Finally,  the  grid  is  shown  in  Figure  4.  If  we  let  the  size  at  point  1  be  0.01.  The 
final  gird  is  shown  in  Figure  5. 


Figure  4.  Grid  of  the  Rectangle  with  uniform  size  required 
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Figure  5.  Grid  of  the  Rectangle  with  non-uniform  size  required 

Example  II 

As  shown  in  Figure  6,  the  L- shaped  domain  is  to  be  meshed.  Four  triangles  are 
used  to  define  the  background.  Two  cases  are  considered.  For  both  cases,  the 
element  sizes  at  nodes  1,  2,  3,  and  4  are  10.  The  sizes  at  node  5  are  10  and  0.1 
respectively.  The  input  data  can  be  given  similarly.  Finally,  the  grids  are  shown  in 
Figures  7  and  8. 


Figure  6.  L-shaped  domain 
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Figure  7.  Grid  of  the  L- shaped  domain  with  uniform  size  required 


Figure  8.  Grid  of  the  L-shaped  domain  with  uniform  size  required 


Example  III 

Finally,  a  rectangle  with  one  hole,  as  shown  in  Figure  9,  is  considered.  The 
boundary  of  the  hole  is  approximated  by  16  line  segments.  Two  cases  for  the 
background  sizes  are  considered  and  listed  in  the  following  two  tables.  The  finally 
grids  are  shown  in  Figures  10  and  1 1. 
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Figure  1 1 .  Grid  for  case  2  of  the  rectangle  with  hole 


Example  IV 

In  this  example,  the  adaptive  femeshing  scheme  provided  in  Section  3  is 
utilized.  The  initial  grid  is  shown  in  Figure  12  which  is  uniform.  After  some  time, 
the  grid  is  regenerated  depending  on  the  error  at  each  element.  During  the 
remeshing,  the  initial  grid  is  looked  as  the  background  grid.  The  refined  grid  is 
shown  in  Figure  13. 


^XX)(XX)Q(XXm)(XX)(XX)(XXX)(J(^ 


Figure  12.  Initial  Grid 


Figure  13.  Refined  Grid 


5.  Conclusions 

The  details  for  the  generation  of  two-dimensional  advancing  front  grid  was 
provided.  The  schemes  for  the  implementation  of  the  generation  were  also 
described.  Two  efficient  seeirching  algorithms  for  locating  the  query  point  have  been 
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presented.  In  the  adaptive  remeshing  the  current  grid  obtained  from  the  advancing 
front  algorithm  is  usually  looked  as  the  background  grid.  Several  examples  were 
provided  to  demonstrate  the  schemes  in  this  report. 
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CONCLUSIONS 

In  the  first  part  P-adaptive  FEM  has  been  developed  and  applied  for  flow  around 
circular  cylinder.  The  second  through  fifth  orders  of  polynomials  are  considered  for 
the  velocity  in  the  adaptive  method.  One  order  less  of  polynomial  is  used  for 
pressure.  Velocity  is  selected  as  the  error  estimator.  A  flow  around  a  circular 
cylinder  with  Re5molds  number  of  1000  was  simulated  using  this  technique.  The 
effects  of  the  highest  order  of  polynomials,  error  tolerance,  and  size  of  the  element 
on  the  accuracy  of  the  drag  and  lift  coefficients  are  surveyed  using  this  flow 
simulation. 

When  the  grid  is  coarse,  the  p-adaptive  technique  is  very  efficient.  It  can  improve 
accuracy  significantly  while  the  computed  time  does  not  increase  much.  For  the 
refined  grid,  the  results  of  the  present  simulation  show  that  the  higher  orders  of 
interpolation  pol5momials,  P=4  and  5  for  example,  do  not  have  much  effect  on  the 
accuracy  of  the  velocity  as  well  as  the  drag  and  lift  coefficients.  The  error  of  the 
velocity  does  not  have  much  effect  on  the  accuracy  of  the  drag  and  lift  force 
coefficients.  If  the  error  estimation  defined  in  this  report  is  reasonable,  the 
aerod5mamic  force  coefficients  are  insensitive  to  the  accuracy  of  the  velocity.  The 
accuracy  of  the  results  improves  a  lot  when  the  prescribed  error  changes  from  10% 
to  5%.  There  is  no  much  change  when  the  error  tolerance  reduces  from  5%  to  2% 
and  then  to  1%  even  though  the  percentage  of  the  elements  with  5th  order 
pol5momial  increases.  The  differences  of  the  aerodynamic  forces  among  the  three 
refined  grids  are  very  minor  while  the  differences  of  the  number  of  equations  are 
very  significant. 

In  the  second  part  h-adaptive  is  investigated  by  first  developing  grid  generation 
program.  The  details  for  the  generation  of  two-dimensional  advancing  front  grid  was 
provided.  The  schemes  for  the  implementation  of  the  generation  were  also 
described.  Two  efficient  searching  algorithms  for  locating  the  query  point  have  been 
presented.  In  the  adaptive  remeshing  the  current  grid  obtained  from  the  advancing 
front  algorithm  is  usually  looked  as  the  background  grid.  Several  examples  were 
provided  to  demonstrate  the  schemes  in  this  report.  Currently  work  is  in  progress 
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to  use  the  grid  generator  to  solve  the  flow  around  circular  cylinder  and  study  the 
issues  in  using  h-adaptive  comparing  to  p-adptive.  Also  work  is  under  progress  to 
develop  robust  solvers.  Some  work  is  published  in  the  following  publication: 

Selvam,  R.P.,  Computational  issues  in  solving  the  incompressible  NS  equations 
(invited  paper),  Proceedings  of  the  2"'^  International  Conference  on  Fluid  Mechanics 
86  Fluid  Power,  P.C.  Jain  et.  al.  (Ed.),  Ajay  Printers  86  Publishers,  IIT  Roorkee,  India, 
Dec.  12-14,  2002  ,  Volume  1,  pp.  1-6 
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